Locating Boosted Kerr and Schwarzschild Apparent Horizons. 
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We describe a finite-difference method for locating apparent horizons and illustrate its capabilities on 
boosted Kerr and Schwarzschild black holes. Our model spacetime is given by the Kerr-Schild metric. 
We apply a Lorentz boost to this spacetime metric and then carry out a 3+1 decomposition. The 
result is a slicing of Kerr/Schwarzschild in which the black hole is propagated and Lorentz contracted. 
We show that our method can locate distorted apparent horizons efficiently and accurately. 



I. INTRODUCTION 

Apparent horizon locators play an integral role in the 
application of black hole excision techniques in the com- 
putational evolution of black hole spacetimes. Excision 
techniques delete the regions of spacetime that contain 
the curvature singularity from the computational do- 
main. Assuming cosmic censorship, these curvature sin- 
gularities are expected to be contained within an event 
horizon. The event horizon is a causal boundary whose 
interior does not causally affect the exterior spacetime; as 
a result it is possible to excise a region within the event 
horizon, thereby excising the black hole's curvature sin- 
gularity. 

In our approach to computationally solving the Ein- 
stein field equations we focus on the use of Cauchy tech- 
niques, in which a 3-f 1 splitting of spacetime into a folia- 
tion of spacelike hypersurfaces, E, parametrized in time, 
is the basis for an evolution in time. The result of this 
splitting is a system of elliptic and hyperbolic partial dif- 
ferential equations in the 3-metric, jij, and extrinsic cur- 
vature, Kij. These are the four constraint equations and 
12 first-order-in-time evolution equations. The Cauchy 
approach starts with an initial spacelike slice with 7^ 
and Kij set by solving an initial value problem (the el- 
liptic constraint equations). One then uses the evolution 
equations to evolve to the next spacelike slice obtaining 
7ij and Kij at the next time (See York 0] for a detailed 
discussion) . 

In the evolution of black hole spacetimes in this manner 
we do not have a complete history of the entire space- 
time and hence do not have a knowledge of the loca- 
tion of the event horizon. Since the event horizon is a 
global object that depends on geometric information for 
all time (or at least until the black hole becomes qui- 
escent) we cannot use it to determine an inner excision 
boundary in our Cauchy evolution. There is, however, 
an alternative, and that is to use the apparent horizon 
surface which is a local object, locatable (if it exists) 



with "fij and Kij at one time. The apparent horizon is 
the outermost marginally trapped surface. It is a closed 
spacelike 2-surface whose future-directed outgoing null 
normals have zero divergence Q|. The apparent horizon 
is slicing dependent and may not necessarily exist even 
though an event horizon does. An example of this is given 
by Wald and Iyer |^ through nonspherically-symmetric 
slicings for the Schwarzschild spacetime. Provided a non- 
pathological slicing is chosen the apparent horizon or any 
trapped surface within it may be used for excising the 
black hole singularity. These surfaces define a local causal 
structure that distinguishes instantaneously escaping null 
rays from those that are certain to collapse. This distinc- 
tion makes their treatment very amenable to computa- 
tional black hole excision techniques. Since these sur- 
faces can be determined with geometric information at 
one instant of time, they are used in practice as an inner 
boundary in Cauchy evolutions. With this purpose in 
mind, we developed a 3D apparent horizon locator that 
utilizes 7ij and Kij on a given spacelike slice of space- 
time and locates an apparent horizon. Once the apparent 
horizons are located, a region contained within the ap- 
parent horizon is excised. Thus the method is really a 
trapped-surface excision. 

There has been a variety of work done on apparent 
horizon location in spherical symmetry, axisymmetry and 
3D. We focus solely on the 3d locators. These can be 
classified into those that use finite difference methods, 
and those that utilize pseudo-spectral schemes. Further, 
one can classify each of these finders in terms of those 
that use fiow methods versus those that directly solve 
the apparent horizon equation either via a minimization 
scheme or Newton's method for root finding. 

One of the first published 3d apparent horizon loca- 
tors was developed by Nakamura, Kojima and Oohara 
Q. Their method expands the apparent horizon shape 
function, r ~ p[6, (p) in spherical harmonics to some max- 
imum I = Irr^aT.- 
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p{e,4>)^Y. E aimYim{9A)- (1) 

i=0 m=-l 

With this expansion Nakamura et al. evaluate the ap- 
parent horizon equation and solve for the coefficients a;„i 
via a "direct" functional iteration scheme. Kemball and 
Bishop Q reimplemented this approach and made mod- 
ifications that led to improved convergence and stability 
behaviour. Anninos et al. |^ and Baumgarte et al. 0] 
implement similar methods that involve an expansion of 
/o(0, 0); the primary differences being that they expand 
in terms of symmetric trace-free tensors and use Powell's 
method for minimization of the square of the apparent 
horizon equation, Eq.(2) below (which is related to the 
expansion of the outgoing null normals). 

Thornburg |^ gives a very good treatise on the use of 
finite differencing to solve the apparent horizon equation 
using spherical coordinates (r, 0, (f) via Newton's method. 
He discusses in general how algebraic Jacobians may be 
applied in a full 3D context. His implementation for hori- 
zon finding however is axisymmetric; his full 3-d finder 
suffers from problems with the z-axis {6 = 0,7r). Our 
method for finding horizons uses closely related concepts 
except that we finite difference in cartesian coordinates, 
eliminating any potential z-axis problems. 

Another class of apparent horizon locators casts the 
elliptic apparent horizon equation into a parabolic one 
as suggested by Tod ||], who suggested the use of flow 
methods in locating apparent horizons. Bernstein px| ] 
implemented Tod's algorithm in axisymmetry using finite 
differences, but encountered problems with differencing 
on a sphere in spherical coordinates in the general case. 

The advantage to flow-methods is that one can start 
with an arbitrary initial guess and flow towards the ap- 
parent horizon(s). In some implementations it is possible 
to find multiple apparent horizon surfaces starting from a 
single initial guess surface (i.e., there is a topology change 
in the course of location of the apparent horizon). Pasch 
PI uses a level-set method to locate multiple apparent 
horizons in 3d. He demonstrates his method utilizing 
time-symmetric conformally flat initial data for multi- 
ple black holes. An hybrid flow/level-set-like method 
utilizing our approach to evaluating the outgoing null 
expansions via Cartesian finite differences has been im- 
plemented by Shoemaker et al. [Q. That method flows 
towards the apparent horizon(s) from an arbitrary initial 
guess allowing for topology changes. Gundlach has 
implemented a "fast flow" method for finding apparent 
horizons. 

In the sections that follow we give a brief discussion of 
the algorithm used and relegate the details to the appen- 
dices. The model spacetime in which all of the results 
are presented is discussed in section III. In section IV we 
discuss tests of the algorithm and demonstrate that for 
distorted apparent horizons for boosted Kerr black holes 
the algorithm fairs well. 



II. BOUNDARY VALUE PROBLEM APPROACH 

On a particular 3D spacelike hypersurface, E, from our 
foliation of spacetime we are given the 3-metric, 7ij and 
the extrinsic curvature, Kij. Let 5 be a closed 2-surface 
in E. At any point p on 5 we can define a spacelike nor- 
mal, s°, to S, and a time-like normal, n", to E. From 
these we can construct the outgoing null normal, /c°, at 
p. If the divergence V ak'^(y a is the covariant derivative 
compatible with the spacetime metric, gab) is zero every- 
where on 5, then 5 is a marginally trapped surface (MTS). 
The apparent horizon is the outermost such MTS. The 
expansion of the outgoing null normals, VqA:" = 0, can 
be rewritten as an equation entirely in terms of quantities 
inE §: 

As* + Kijs'^s^ -K = 0. (2) 

Here Di is the covariant derivative compatible with the 
3-metric jij and K is the trace of the extrinsic curvature. 

The apparent horizon equation is an elliptic partial 
differential equation on S (for a function of coordinates 
on S). This can be made apparent by noting that an 
MTS is a closed 2-surface; spherical coordinates are a 
natural set of coordinates for S. The location of S can 
then be written as the radial distance from the origin of 
the coordinate system, r = p{9,(f>). In general one can 
generate a foliation of such closed spacelike 2-surfaces 
based on the distance from the MTS. This is given by 

ip = r-p{9,cP), (3) 

where the (p — level surface is the MTS. From (p we 
define a spacelike vector field, the normal 

= Y^djp/y^j''''dk(pdip, (4) 

at every point on these level surfaces. Substituting this 
into Eq.(^ results in a second order elliptic partial dif- 
ferential equation on S, 

F[p] - r'dadbV + r\db^ ~ \u:-^r'^l^dbpd,^da^ 

+Lo-iKabr'-f'"'d,pddip-Lo^K ^0, (5) 

where u = ^'^'^dc^pddf and P^^^ is the connection coeffi- 
cient associated with the 3-metric ^ab- 

Our approach involves casting Eq.(^ as a bound- 
ary value problem on S. As stated, points on S are 
parametrized in spherical coordinates {9 £ [0, 7r],(/) £ 
[0, 27r)). S is discretized into an uniform mesh, 5, of 
No X Nij, points where Ng — — Ng. The domain 
on 5 is (0 < 6* < tt; < (/) < 27r) where at the poles, 
6 = O. TT all points are identified as one. The (j> = 2tt 
branch cut is identified with the = line. The bound- 
ary conditions simply are periodicity at ^ = 27r and (j) 
identification at 6' = 0, tt. These boundary conditions 
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are key to avoiding the coordinate singularities at the 
poles in combination with using Cartesian coordinates 
to discretize partial derivatives on S. We treat as a 
function of Cartesian coordinates x, y, z and center on 
each mesh point of 5 a 3-d Cartesian difference stencil of 
27 points. Using the form Eq.(|^) we interpolate values of 
Lp{x^y, z) onto each of the 26 stencil points surrounding 
each S stencil point. (See the appendix for more details.) 
Using this difference stencil we can evaluate first, second 
and mixed derivatives of 1^9(2;, y, z) as required by the dis- 
cretized version of Eq. . At every point on S we then 
construct the residual F[p] on S (Note that the discrete 
version of a continuum quantity, T, is denoted by, T). 

The problem at hand is to solve for a p that yields 
F[p\ = 0. Since F[p] is a nonlinear operator (as shown in 
Eq.(||)), we use Newton's method to solve for p. Given 
an initial guess surface, p = po, we wish to find a 5p (the 
change in the surface) that leads to F[po + dp] = or, to 
lowest order. 



F[po + Sp] - F[po] 



dF[p] 



dp 



dp + O{Sp^)=0. (6) 



P=PO 



The Jacobian of F[p] is defined to be 

dF 



J 



dp 



(7) 



In the discretized case, j is an x iV matrix, where N 
is the total number of points used in the discretization. 
To obtain a Sp that leads to F[p + Sp] = we must solve. 



J-6p = -F[p] 



(8) 



for Sp. 

Computationally our tasks are to first evaluate the dis- 
crete form of the Jacobian matrix, J and second to solve 
the discrete form of Eq.(||). We numerically compute 
the Jacobian matrix by perturbing the surface pointwise 
and examining the effect of of the perturbation on the 
residual, F. Let p denote "independent" points in the 
computational mesh, S. By independent we mean the 
unique points (points modulo boundary identification) 
on S. In particular from the identifications made ear- 
lier, there are N'^ — 2Ns + 2 independent points in S. 
Ng = Ng = points at each of the poles are treated 
as one point, ji = 1 represents the 9 = point for all 
the points (0 < (j) < 2tt) and p, = - 2Ns + 2 is 
the = n point. Eq.(|8|) then becomes a linear system of 
equations where J is a {N^ - 2Ns + 2)x {N^ - 2Ns + 2) 
matrix and F and Sp are vectors of length — 2Ns + 2. 
The pD component of the Jacobian is then computed by 
perturbing p at the v-th point and computing the change 
in the residual, F at the p-th point. Using a first order 
forward difference approximation we have. 



where e is the amount by which we perturb the surface. 
We define e to be the perturbation parameter. The pro- 
cess for generating the components then involves numer- 
ically evaluating F in only a small neighborhood of the 
h'-th point since F has a domain of dependence depen- 
dent on the finite difference operators used. In this case 
the operators are (finite difference) derivative operators 
convolved with interpolation operators. 

The solution of Eq.(||) is achieved via Newton's 
method. We solve for the change Sp that leads to a new 
surface r = p^.{9,(j)) that yields F[p^] ^ up to 0{h?), 
where h is the Cartesian stencil spacing, proportional to 
SO. Newton's method then involves updating p[6, (p) by 
Sp. li F were a linear operator then one iteration would 
result in a Sp that leads to a solution. Since F is nonlin- 
ear we have to iterate until the _L2-norm of the residual, 
II-FII2, is driven down to a chosen stopping criterion. We 
discuss the implementation details in the appendices and 
discuss further the properties of the Jacobian and New- 
ton's method in the results section. We now turn our 
attention to the model spacetime in which we shall con- 
duct our numerical experiments. 



III. 3+1 SPLITTING OF THE KERR-SCHILD 
METRIC 

In the rest of the paper we focus on tests of the al- 
gorithm based on boosted Schwarzschild and Kerr black 
holes. The particular form that we use is given by the 
Kerr-Schild metric: 



9tiu = Vnv + '^Hl^l^ 



(10) 



(9) 



where is an ingoing null vector (i.e: g^^lf^lu ~ 
rj^'^lfjli, = 0), H is a. scalar function of the spacetime 
coordinates and 77^^ is the Minkowski spacetime met- 
ric. This metric describes the Kerr and Schwarzschild 
spacetimes. We note that under a Lorentz transforma- 
tion the spacetime metric is form invariant. By definition 
such a transformation takes 77^^ — *■ r]^^, and and H 
are transformed to a new null vector and left unchanged 
(though evaluated at the new coordinate labels for the 
same event) respectively. This property makes our anal- 
ysis easier since a 3-1-1 decomposition of Eq.(|l^) has the 
same form as a 3-f 1 decomposition of the boosted metric. 
As we shall see, we only need to specify H, and their 
spacetime derivatives in order to obtain the 3-metric and 
extrinsic curvature on S. 

For a vacuum spacetime, is geodesic and in the Kerr 
and Schwarzschild spacetimes is the tangent to geodesies 
of ingoing photons. The null nature of Eq.(|l^) leads to 
a slicing of these spacetimes that is well behaved at the 
horizon. That is, spacelike slices penetrate the horizon 
and hit the black hole singularity. This is a desirable 
property for black hole excision in computational appli- 
cations and this metric has shown itself to be a good 
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choice for the study of single and muhiple black hole 
evolutions with exicision. 

For the Kerr spacetime, H and are given by 

Mr^ , , 

(11) 



H 



and 



rx + ay ry ~ ax z 



where r is given by 



x"^ + 



(12) 



(13) 



or 



= i (p2 - a^) + (p2 _ „2)2 + „2^2. (14) 

M is the mass of the Kerr black hole and a = J /M 
is the angular momentum of the black hole and p — 

In the a — > limit we get the Schwarzschild metric in 
ingoing Eddington-Finkelstein coordinates where 

M , , 

H=—, (15) 



(16) 



and r — ^/ x"^ + y'^ + z^ . 

In a spacelike slice of either Kerr or Schwarzschild 
spacetimes the apparent horizon is known to coincide 
with the intersection of the event horizon with that slice. 
In the Kerr spacetime then the apparent horizon is a sur- 
face of radius r = r i : 



and area 



A = 47r(r?_ + a^ 



(17) 



(18) 



In the more general nonstationary case the apparent hori- 
zon and event horizon will not coincide. We use the prop- 
erties of the Kerr and Schwarzschild spacetimes to test 
out our method for finding horizons. 

To get the spacetime metric for a boosted black hole 
consider O to be the rest frame of the black hole, with 
coordinates (i, x'). Let O be another stationary frame 
with coordinates {t, x^) such that O is related to O via 
a Lorentz boost along the v — {vx,Vy^Vz) direction: in 
the O frame the black hole moves in the v direction with 
boost velocity, v [Sijffv^ = 1). As usual, we define 7 = 
1/Vl — v"^ . H(xji) and Ifi (bar denoting O frame) now 
transform as 



Hix^) = H{Klx,) 



and 



l^=kll-,{KZ^x^). 
These preserve the form of (p^. 



(19) 



(20) 



A. 3+1 Decomposition 

The standard ADM 3-1-1 form of the spacetime metric 
is given by 

ds^ = ^a^dt^ + {dx' + P'dt) {dx^ + 13' dt) (21) 



If we compare (|1C|) to (|21j) and use the property that 
l^lf^i = 0, we find that the lapse is given by 



1 

v/1 + 2Hl^ ' 
and the shift is given by 

Pi — ^Hlth 

or 

= 2HltS'Hj/{l + 2Hlj). 
The 3-metric is given by 

lij = Vi] + 2Hlilj. 



(22) 



(23) 



(24) 



(25) 



as expected and the extrinsic curvature is determined 
from 

K^J = -^a^J/2a + D,(3j + Dj(3, (26) 
- -dt{Hklj)/a + 2 {D,{Hltlj) + Dj{Hltli)) (27) 



and 

^ s^J _ 2H5'^5i^lilk/{l + 2Hl1). 

Note that 



(28) 



det'-i.j = 1 + 2Hl1. (29) 

To obtain the 3-metric and extrinsic curvature we need 
to specify H , lf^,duH and d^li, and substitute these into 
Eq.(p5|) and Eg. (p7|) . In order to evaluate df^li, and df^H a 
specific choice of spacetime has to be made. For example 
for the Kerr spacetime we take the expressions for H and 

from Eq.([l2|) and Eq.(|ll|), compute their derivatives 
and substitute. This gives us a "Kerr-Schild" slice of the 
Kerr spacetime. 



IV. RESULTS 

In the presentation that follows we first conduct a se- 
ries of basic tests of the algorithm using metric and ex- 
trinsic curvature data just presented. In the second part 
of this section we set up a 3-dimensional Cartesian grid, 
S, of points, on which we define a coordinate sys- 
tem where the black hole (either Kerr or Schwarzschild) 
is placed at the origin. Using Eq.( p5| , p7| ) again we gen- 
erate "fij and Kij on S everywhere but the region that 
contains the curvature singularity (for Schwarzschild at 
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and Kerr p = 



z2 < 



With the data on S, we start the apparent horizon lo- 
cator with an initial guess surface which is a 2-sphere of 
radius rg. The horizon locator surface mesh sizes used 
in this calculation are Ns = 33,41,49,65,81,97. The 
stopping criterion for the Newton iterations was deter- 
mined empirically to be 10^^. The cases we present are 
(1) ti = 0, a = (unboosted Schwarzschild), (2) w = 0, 
a ^ (unboosted Kerr), (3) w 7^ 0, a = (boosted 
Schwarzschild) and (4) f ^ 0, a 7^ (boosted Kerr). 



A. Tests with Eddington-Finkelstein metric data 

First, with = 0, Uj, = 0, 1)2 = and a = (un- 
boosted black hole in ingoing Eddington-Finkelstein co- 
ordinates) we show some basic tests of the apparent hori- 
zon locator. Most importantly we show that solutions 
obtained with our locator are 0{h?\ With these data 
all components of 7^ and Kij are non-zero. The lat- 
ter property makes this a good initial model problem to 
work with, because the computation is fully exercised in 
an analytically tractable situation. As stated earlier the 
apparent horizon is expected to be located at r = 2Af . 
These tests are conducted with data specified analytically 
where required. 



1. Residual Evaluation and Second Order Convergence 

We place the black hole at the origin of the computa- 
tional domain(a: = 0, y = 0, 2: = 0). In spherical symme- 
try for this metric the apparent horizon equation becomes 
the algebraic equation. 



1 - 2M r 
F r = = 0. 

ry/l + 2M/r 



(30) 



A plot of F{r) is shown in FIG.(|lj). At r = 2M we have 
F = 0. A useful test of the evaluation of the expansion of 
the outgoing normals F{r) is to see if indeed the residual 
F[p] is correctly evaluated to 0{h?) as 



F = F + 63/1^ 



(31) 



where 62/1^ is the leading order truncation error term. 
Given that the exact value is known for F[p\ we can 
approximate the leading order truncation error. We 
carry out a convergence test by evaluating F[p\ on a 
2-sphere of r = 2M for a series of mesh sizes, Ng — 
17,25,33,49,65,96,129. We examined the behavior of 
log ||F||2 (where this is the L2 norm) versus logiV^, where 
Ns is the number of mesh points on one side of the 
Nf, X TVs mesh. At r = 2M F ^ 62/1^ and so the L2- 
norm, ||F||2 ~ ||e2|j2/i^- Since h oc l/Ng we expect that 
if the residual is 0{h?) then the slope of a plot of log ||F||2 
versus logA^s should be —2.0, which we validated via a 



least squares fit. A closely related test is to also evalu- 
ate log Hp — p\\2 versus logiV^, where p is the numerical 
solution from the apparent horizon locator and p is the 
exact horizon location. FIG.(^) shows the result. From 
a least squares fit to a straight line the slope is found to 
be —2.1 which validates our solution as 0{h^). 



2. Jacobian 

For the same 2D mesh discussed we generate the Jaco- 
bian matrix for a single Newton step. FIG.(^ shows the 
structure of the matrix for a 33 x 33 run. There are 1025 
independent points on S and hence Jpp is a 1025 x 1025 
matrix. The dots in the figure are non-zero Jacobian 
entries. There are seven bands in this matrix with 2 ad- 
ditional ones in the vicinity of the poles at fl — 1 and 
p, = 1025. The structure reflects the domain of depen- 
dence of the finite difference operators used in the eval- 
uation of F. Here it comes from a combination of inter- 
polations and Cartesian finite differencing. Near p = 1 
and p — 1025 the additional bands come from our special 
choice of interpolation stencils at the poles, as discussed 
in the appendix. 

The structure reflects the fact that a perturbation at a 
single mesh point affects the residual in a small neighbor- 
hood around it so we can optimize the generation of the 
Jacobian to 0{N) by evaluating F[p + Sp] only in a small 
neighborhood of the perturbed point. The Jacobian gen- 
eration was found to be order 0{NP) where p — 1.08 and 
N is the total number of independent points on the 2D 
computational mesh. 

A matrix A is defined to be diagonally dominant p^ ] 
if its elements, Aij, satisfy 



^\A,j\<\ Au I for all i. 



(32) 



3=1 



We found that the Jacobian is not diagonally dominant 
since the inequality in Eq. (|3^ ) is not satisfied for all i and 

This is of interest since for some iterative solution tech- 
niques (Gauss-Seidel and SOR for example) a sufficient 
condition for the solution of a linear system, A - x — b, is 
that the matrix, A, be diagonally dominant. In our case 
we concluded from early experiments that indeed such 
simple iterative solvers did not converge for this prob- 
lem. 

The Jacobian matrix is not symmetric but it is well- 
conditioned for the spacetimes that we have considered. 
For a 33 X 33 run the Jacobian has a condition number 
of K of about lO'' to 10^ where. 



l^llll^- 



(33) 



(In their definition of the condition number, Dongarra et 
al. |l6| use the Li-norm.) The condition number tells us 
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how close the matrix A is to being singular. A very large 
condition number or a reciprocal condition number close 
to machine epsilon tells us that A is singular. An identity 
matrix has a condition number of 1. To estimate k we 
used the LINPACK library routine, DGECO. 

3. Solution of the Linear system 

As stated before to locate the apparent horizon using 
our technique we have to obtain a solution, 6p, to the 
linear system 

J-Sp= -F[p] 

which is the discrete form of Eq.(^). 

Since the properties of the matrix do not allow us to use 
the standard iterative methods such as Jacobi and Gauss- 
Seidel methods. We use a modified conjugate-gradient 
method due to Kershaw (The standard form of the 
Conjugate gradient method will not work since J is not 
symmetric.) Kershaw's method, termed the Incomplete 
LU-Conjugate gradient method (ILUCG), can solve any 
linear system, A ■ x = b, with A being any nonsingu- 
lar, sparse matrix. The method involves preconditioning 
the matrix via an incomplete LU decomposition. This 
method has worked quite well for our purposes. In prin- 
ciple other schemes for solving the resulting system can 
be used. 



4- Solution for the apparent horizon location 

Using the ILUCG method to solve for Sp we demon- 
strate the Newton solver's ability to sucessfuUy locate 
apparent horizons in the Eddington-Finkelstein metric 
data. The apparent horizon in this data is a 2-sphere of 
radius 2M. Using a 2-sphere of radius ro, centered on the 
origin as the initial surface we carried out a series of runs 
for ro = 0.5..3.0 with a 33 x 33 mesh. Table (||) shows the 
radius, rp, of the initial starting 2-sphere, the number of 
Newton iterations taken if it converged, the final resid- 
ual value and the solution error. The stopping criterion 
is that the norm of the change in the solution, \\Sp\\2 be 
less than 10~^°. We see that for all the cases, provided 
the solver managed to drive \\Sp\\2 below 10~^° the final 
percentage error remains fixed; what differed in each case 
was the number of iterations taken and rate of conver- 
gence. Once the solver has driven p into the vicinity of 
the solution the Newton convergence is quadratic. This 
happened in the 6th iteration in the rp = 0.5 run. For 
the To = 2.0 run the solver took four iterations to con- 
verge down to the stopping criterion. In another series 
of runs with Ng = 33, M = 1 and ro = 2.5 the pertur- 
bation parameter, e, was varied from 10~^ to 10~^. An 
optimum value of e for this metric data was found to be 
about 10~^ to 10~^. In general the stopping criterion 
need not be as stringent as we have set it. In numerical 



TABLE I. This table shows a series of runs carried out for 
various initial conditions. The initial surface is a 2-sphere of 
radius ro where ro is shown in the first column. The number 
of iterations taken for the solution to achieve ||5p||2 < 10"^" 
is given. Note that the final error in the solution remains 
a constant provided the solver is able to drive ||(5p||2 below 
the specified stopping criterion. The perturbation parameter 
used to generate the Jacobian was 10~^. 
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3.00 


12 


9.9E-11 


3.0E-Q5 


1.5E-Q3 



spacetimes where the metric data will have truncation 
error associated with them the truncation error of F is 
expected to be much larger than our test stopping cri- 
terion. In that case a larger stopping criterion should 
to be chosen to avoid wasting computational effort. The 
optimum value is dependent on the error in the metric 
data. In the results we present in the paper however we 
drive the residuals down as far as possible. On the other 
hand, a perturbation parameter e must be chosen such 
that F[/5-|-e] — is sufficiently large that this expression 
is not dominated by truncation error. 

5. Numerical Metric data 

Tests described so far used data analytically computed 
at each point as needed. Since the ultimate goal is to in- 
corporate this apparent horizon location algorithm into 
an evolution code it is useful to gauge the performance of 
the algorithm with numerical metric data and with the 
data structures expected in the real application, where, 
for instance, part of the domain is excised from consid- 
eration. Thus we set up the same Eddington-Finkelstein 
data on a 3D Cartesian grid of points, with a re- 
gion of this grid excised to emulate the situation in an 
evolution code where the interior of the black hole is ex- 
cluded. The apparent horizon surface which is embedded 
in this 3D Cartesian grid typically does not lie on Carte- 
sian grid points and as a result an interpolation tool is 
required. If the surface mesh, during the course of the 
Newton iterations, overlaps the excised region then ex- 
trapolation is required. We make use of an interpola- 
tor/extrapolator written by S. Klasky |^ to obtain the 
3-metric, extrinsic curvature and the spatial derivatives 
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of the 3-metric at any point. The use of an interpolator 
brings in truncation errors associated with the interpola- 
tion/extrapolation operations. In the following we show 
that even with extrapolation errors the solver works quite 
well in locating apparent horizon surfaces. 

We set up an uniform 3D Cartesian grid, S, of size n?" . 
On this grid we excise a region interior to a sphere of ra- 
dius Rm centered at (xm, J/m, ^m) so that the metric data 
is defined for r > Rm and undefined for r < Rm where r 
is the Cartesian distance in Kerr-Schild coordinates from 
the excision center. 

In the following discussion on radial and offset appar- 
ent horizon locations we take n = 65 for the Carte- 
sian grid (with h = 1/8) and Ng = 33 for the surface 
mesh. The stopping criterion used in the horizon finder 
is (3 = 10"^. That is, if \\Sp\\2 < /? then the Newton it- 
erations are stopped. The perturbation parameter, e, is 
taken to be 10""*. The interpolator tool is used to fourth 
order. The initial guess surface used is a sphere of radius 
ro — 2.1M centered at the origin of the Cartesian grid. 
With these parameters we carry out two set of tests. The 
first is a radial test of the horizon locator with the use of 
the interpolator and the second is an offset test. These 
tests examine the effect of extrapolation of metric data 
on the residual, F, and solution to the apparent horizon 
equation. 

6. Apparent horizon location (Radial tests with excision) 

In the first test case we center the black hole at (0, 0, 0). 
The masked region is also centered at (0, 0, 0). We carry 
out a series of tests with the excision radius, Rm, vary- 
ing from 1.5 to 2.6. Thus the apparent horizon is in the 
defined region {Rm < 2M) for some of the tests, and for 
others it is inside the excised region {Rm > 2M). This 
provides evidence of the effect of extrapolations on the 
residual of the apparent horizon equation, F, and the er- 
ror in its solution. FIG.(Q) illustrates the behaviour of 
the L2-norm of the residual, ||-F||2, as a function of Rm- 
FIG.Q shows the percentage relative error of the solu- 
tion of the apparent horizon equation as a function of 
Rm ■ The percentage relative error is calculated from the 
exact solution for r — p{9, cj)) = p = 2M. With this ex- 
act solution, we calculate the percentage relative error, 
e EE Hp - p\\/p X 100%. For Rm < 2M the derivative 
interpolator/extrapolator uses interpolation for regions 
near the apparent horizon location (?' = 2M), while for 
Rm > 2M it uses extrapolation. As Rm increases fur- 
ther the errors due to extrapolation increase, as expected. 
This can be seen in FIG.(|) where \\F\\2 increases quickly 
for Rm ^ 2.4, as does the error shown in FIG.(|). At 
Rm — 2.5M, the solver could not bring ||i5p||2 down to 
below 10^**, and so failed to meet the stopping criterion. 
This can be understood in terms of the Cauchy-Schwarz 
inequality |jl^. Since || J • Sp\\ ~ we have that 



II^PII > (34) 

At Rm = 2.2 where ||F|| ~ lO'^ and ~ 10-^ we 
have from Eq.(|3|) that |1 J|| - lO^. Therefore at Rm = 
2.5M we expect with ||F|| - 10"^ that \8p\\ - lO"'^. 
By relaxing the criterion past Rm ~ 2.5 we can still ob- 
tain a solution. Past Rm — 2.6 the convergence progres- 
sively worsens. For example, at Rm = 2.9, ||i^|| could 
not be brought below 10"'^, and the solution error is 
5%. The amount of error sustained from interpolation 
of the metric data is dependent on the resolution of the 
Cartesian grid and the behaviour of the functions being 
interpolated. If the gradients of 7^ and Kij are large 
near the horizon then a larger interpolation error is sus- 
tained. This in turn leads to a larger truncation error in 
F . In the numerical evolution of black hole spacetimes 
with excision then buffer zones may not be necessary for 
the location of apparent horizons. However, for other 
reasons buffer zones might be necessary. 

7. Locating Offset apparent horizons 

We examine behaviour of the locator with the deriva- 
tive interpolator for a black hole offset so that it overlaps 
the excised region. This is important in tracking moving 
black holes [pO| . 

The center of the masked region is at (0, 0, 0) 
and the black hole of radius 2M is centered at 
{S/VS, S/^/S, S/VS), so that the radial distance between 
the mask center and the hole is S. With a grid spacing of 
h = 1/8, an offset of 5 = 1 corresponds to approximately 
8 grid zones. FIG.@ shows the percentage relative er- 
ror in the apparent horizon location as a function of the 
offset S. As the graph illustrates, up to (5 = 0.7 the per- 
centage relative error is below one percent. (At 6 = 0.7 
the percentage error is 0.6%.) From S — 0.7 onwards, 
however the solver becomes sensitive to initial conditions 
and extrapolation errors and quickly ceases to converge. 

At 5 = 0.7, about 5-6 grid points offset, we are still 
able to find horizons. Generally in explicit time-evolution 
codes the CFL condition |^ restricts the black hole mo- 
tion from one time slice to another, to be less than one 
zone {S < h or about 5 ~ 0.1 in our test case). Hence we 
expect, based on the results for our model spacetime as 
shown in FIG.(^, that in such an evolutionary scheme 
with a similar resolution we will be able to track black 
hole apparent horizons to less than 0.1%. 

B. Apparent Horizons in Boosted Kerr data 

In this section we now focus on apparent horizon loca- 
tion for boosted Kerr and Schwarzschild black holes. For 
the data that follow we excise a 2-sphere of radius r > a 
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centered about the origin from the computational do- 
main to avoid the ring singularity structure of the Kerr 
black hole. Using the interpolator tool in conjunction 
with the apparent horizon locator, we locate horizons for 
various values of the angular momentum parameter, a. 
The horizon locator begins with a trial surface which is 
a 2-sphere of radius rg = 2M. The locator was run for 
a = 0.0, 0.1, 0.2,..., 0.9 at Ns = 33. FIG.(0) shows a 
cross-section of the horizon in the xz plane as a function 
of a. The apparent horizon is seen to have the shape of 
an oblate spheroid. In the runs used to generate these 
data we used e — 10~^ and a stopping criterion that en- 
sured that the ^2-norm of F on the computational mesh 
was less than 10^^^. 

FIG.(H) shows the Z2-norm of the error in the solu- 
tion, ||(f — r+)/r+||2, versus mesh size. This set of runs 
was carried out with a = 0.9 and Ng = 17, 25, 33, 49, 65. 
Where r+ is given by ( |l7| ) and the f is computed from 
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(-2 2\2 , 2 2 

\p — a ) + a z 



1/2 ■ 



1/2 



(35) 



where p — i/a;^ + + is the solution from the appar- 
ent horizon locator. From a least squares fit the slope is 
found to be —2.1 and the solution is 0{h?). 

The area of the event horizon in the Kerr spacetime 
p3| is given by 



4:7: (r 



(36) 



Let A be the computed apparent horizon area. FIG.(||) 
shows the percentage errors ||(^ — A)/ A x 100%||2 ver- 
sus a for various resolutions Ns — 17, ... , 65. The area of 
the apparent horizon is computed via a technique which 
projects the 3-metric, jij, onto the 2-surface to obtain an 
area element \/(^^dOd(j), and then computes the surface 
integral. FIG.(^ shows the percentage errors in the area 
for increasing resolution. We now consider Schwarzschild 
and Kerr black holes boosted in the yz-direction. That is, 
we look at Vx — 0,Vy — 1/^/2, Vz = l/\/2 and a = 0,0.9. 
In each of these cases we locate apparent horizons for 
V — 0, 0.1, . . . , 0.9. From v = Otov = 0.8 we started with 
a two-sphere of radius 2M and found an apparent hori- 
zon with outgoing expansions driven down to 10~^^. For 
w > 0.8 we had difficulty driving the expansions down. 
As a result we utilized the solution at v = 0.8 as an ini- 
tial guess and were subsequently able to find horizons by 
stepping every 0.25 from v — 0.8 to t; = 0.9. We used 
e = 10^^ again for these runs. At u = the initial guess 
is the apparent horizon and there within the six New- 
ton iterations the expansions were driven down around 
10~^^. The first Newton iteration took the expansions 
down around 10~^. For v = 0.5 starting from an initial 
guess of a sphere of radius p = 2M it took four Newton 
iterations to drive the expansions down around 10"^ and 
nine Newton iterations to get down to 10~^^. Typically in 
a numerical time-evolution of such a spacetime we would 



not need to drive the expansions down to this level. If we 
are utilizing a surface within the apparent horizon as an 
excision boundary then we need only to drive the expan- 
sions down far enough to be certain an apparent horizon 
is present. FIG. (00) shows the yz-cross-section of the 
apparent horizon for various boost velocities compared 
against an unboosted black hole apparent horizon cross- 
section. We find that the apparent horizon is Lorentz 
contracted in the yz-direction in the boosted coordinates. 
We have considered a slice of such a boosted spacetime 
in which the event horizon appears Lorentz contracted in 
the resulting coordinates. We know that in these space- 
times the apparent horizon should coincide with the event 
horizon and we find that this is indeed the case. First, 
the area of the apparent horizon coincides with the area 
of the event horizon which is invariant under a boost. 
FIG.(pTl) shows the error in the apparent horizon area 
as a function of v for various resolutions. Wc find that 
with increasing resolution the error in the area converges 
towards zero. This demonstrates that the area of the ap- 
parent horizon is invariant under a Lorentz boost. This 
is coupled to an interesting property of the Kerr-Schild 
type of metrics that r = r+ remains fixed for Kerr and 
Schwarzschild black holes. This is illustrated in FIG.(p^) 
where we show the error in the radial coordinate r = 2M 
on the apparent horizon for various boost velocities. In 
this case the black hole is boosted in the a;yz-direction for 
generality. That is, = l/VS,Vy = l/VS, = l/VS 
and a = 0. Here r is computed from the boosted coor- 
dinates. We find that r converges towards 2M for in- 
creasing resolution satisfying yet another property of the 
boosted Kerr-Schild spacetime. 

FIG. (13) shows surface plots of the apparent horizon 
for V = 0, 0.3, 0.6 and v = 0.9 displayed in Kerr-Schild 
Cartesian coordinates. Note how distorted the apparent 
horizon gets with increasing boost velocities. As seen in 
the figures for the j/z-boosts the boosted apparent hori- 
zons in this case are always contained within the apparent 
horizon for u = 0. That is, the boost contracts the ap- 
parent horizon in the boost-direction. Again for a boost 
velocity oi v — 0.5 it took the solver eight Newton steps 
to drive the expansions down around 10~^^. The stop- 
ping criterion used in this run was 10~^^ and the final 
expansions are ~ 10"^"^. On average it took four Newton 
steps to drive the expansions down to ~ 10~^ and nine 
Newton steps to 10^"'^^ starting from an expansion of 
0.1. 

In the case of a boosted Kerr black hole with a — 
0.9 the results are again very similar to those of the 
Schwarzschild black hole. Note that now with a = 0.9 
and u — > 0.9 we get even more distorted apparent hori- 
zons. These results show that this algorithm for find- 
ing apparent horizons does quite well with such large 
distortions. In addition the cost of finding these sur- 
faces increases by only two additional Newton steps. 
FIG.(p^ shows the yz-cross-sections for the apparent 
horizon found for a = 0.9 as a function of v. Again 
the boosted apparent horizon is contained within the un- 
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boosted one and Lorentz contracted. FIG.([l5|) shows the 
error in the area for the same data. With a = 0.9 we ex- 
pect that the area should be '--^ 36. The graph shows that 
for increasing resolution the error in the area tends to- 
wards zero. Hence the area remains fixed with increased 
boost velocity as is expected. 

Similarly, r computed from the boosted coordinates re- 
mains fixed at r+ as is shown by FIG.(^6|). The apparent 
horizons found here were obtained with Vx — Vy = 

l/\/3, Uz = l/VS and a = 0.9. That is, the boost was 
in the xyz-direction with magnitude v. Again we find 
that r on the apparent horizon converges to 1.4 
with increasing resolution for all boost velocities. At 
a resolution of 33 x 33 we have an percentage error of 
8% and 1% at 65 x 65. FIG.® shows surface plots 
of the apparent horizon for the boosted Kerr black hole 
for V — 0,0.3,0.6 and 0.9. Note how distorted the fi- 
nal apparent horizon surface is. Our algorithm required 
one more Newton step to drive the expansions down to 
10"^"^ for t; = compared to v = 0.9. It took six Newton 
iterations to drive the resolution from about 0.2 at the 
inital step to 10~^ for both boost velocities. Hence, this 
algorithm has the advantage that given sufficient resolu- 
tion on the computational mesh, the work done does not 
drastically increase for increasing distortions. 

V. DISCUSSION 

We have demonstrated in this paper that our method 
based on finite difference techniques is viable for locating 
very distorted boosted Kerr black hole apparent hori- 
zons. We have shown that the located horizons obey 
the expected analytical rule, of invariance of the area of 
the event horizon, in cases corresponding to at-rest or 
boosted single black holes, where the apparent horizon is 
known to coincide with the event horizon. We have ad- 
ditionally given a number of computational tests demon- 
strating the behavior of the tracker on interpolated or 
extrapolated data which is realistically similar to that 
from evolutions. In other contexts algorithm has been 
thoroughly tested with the canonical set of test problems 
such as the two and three black hole initial data sets and 
additionally in an evolution code tracking the apparent 
horizon for a Schwarzschild black hole in geodesic slicing 
pl| and demonstrated to be capable of tracking apparent 
horizons in boosted Schwarzschild data . Those tests 
and the tests given here show its viability as a method for 
locating black hole apparent horizons and using them for 
black hole excision. Since black hole excision is essential 
for long-term evolutions of single or multiple black hole 
spacetimes. It is very useful to have efficient apparent 
horizon locators that can locate apparent horizons "fast" 
relative to the time taken for an evolution time step. 

Our algorithm is dominated primarily by computations 
of the Jacobian matrix in the use of Newton's method. 
These operations are optimized such that they scale as 



0{N) approximately where A'' is the total number of 
points on the two-dimensional mesh used for the solution. 
With a 33 X 33 mesh we find that each Newton iteration 
takes on average 20 Origin 2000 CPU seconds. This is 
independent of the distortion of the apparent horizon. 
However the number of Newton iteration steps is deter- 
mined by the "distance" of the initial starting surface 
from the final solution. During the course of an evolu- 
tion it is expected that the apparent horizons over several 
timeslices will be "close" enough to each other that two 
to three Newton iterations will be sufficient to locate the 
horizon at low accuracy with the expansion of the outgo- 
ing null rays on its surface being at the level of lO^'^ or 
10^^. Obtaining a better accuracy requires more New- 
ton iterations and the number of timesteps taken depends 
also on the accuracy of the background metric data. Typ- 
ically in our model problems eight iterations will take us 
below 10-1°. 

One of the drawbacks of a Newton's method for find- 
ing apparent horizons is its sensitivity to the initial guess. 
An initial guess outside of the radius of convergence will 
not lead to a solution. Additionally Newton methods 
are known to be sensitive to high frequency components 
in the solution. This is demonstrated in axisymmetry 
by Thornburg Sensitivity to the initial guess can be 
easily handled by combining the Newton method algo- 
rithm with apparent horizon trackers that are based on 
flow methods. The flow finder is used to obtain an initial 
guess for the Newton method which then converges on 
the solution very quickly. 

The efficiency of our boundary-value method can be 
compared to the efficiency of other approaches (variations 
of flow mothods) due to Tod, as developed by Shoemaker 
et al. ||l^ and fast flow methods developed by Gundlach 
p3[ . The flow method is based on a parabolic partial 
differential equation whose rate of convergence to the so- 
lution slows as it approaches it. Typically for a 33 x 33 
run the flow method takes on the order of thousands of 
seconds to converge down to expansions of 10^^. The ad- 
vantage of this method however is its ability to flnd mul- 
tiple apparent horizons from an arbitrary initial guess. 
Combined with the Newton finder this will result in a 
robust apparent horizon finding scheme. 

We can also compare the effort to spectral decompo- 
sition methods. We concentrate on a method similar to 
that of Nakamura et al. Q in which the equation for the 
apparent horizon surface is written: 

p{e,q^) '^ImYlmiO,^). (37) 

1=0 m=-l 

We do not have access to an apparent horizon finder 
based on pseudo-spectral methods but we will analyti- 
cally compute the coefficients for the case of a boosted 
Schwarschild black hole; this will give some insight into 
the range of harmonics required, and some idea of what 
scaling of these methods might be. 
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In Kerr Schild coordinates, the hole, with boost in the 
z-direction, has the shape of a spheroid. 



where 



z 



1, 



(38) 



where we have supressed the y-direction. 

Notice that the axes a,b of the eUipsoid obey &^/a^ — 
1 — w^, which demonstrates that the eccentricity is di- 
rectly proportional to the boost velocity, e — v, for this 
case. Hence even ellipses with moderate ratio of axes, 
such as that for v=0.9, where the ratio is a little less 
than 0.5, have moderately large eccentricities. We will 
approximate the form Eq.(pq) with an axisymmetric se- 
ries of the form Eq.(^) (the general case would have 
nonaxisymmetric terms also). We find it more conve- 
nient to work with Legendre polynomials than with the 
spherical harmonics directly. 

Since we work with Legendre polynomials, we drop the 
y- coordinate in the spheroid expression, to obtain : 



R{9) =&/vl-e2sin^ 



= 5/Vl-e2(l-g2) 



(39) 
(40) 



where q = cos 9. 

To obtain the expansion of expression Eq. ( |40| ) in terms 
of -Pm? we first expand using the binomial theorem. 



OO / /n \ 



(41) 



This converges for all w < 1. 

Using the binomial theorem again for — we sub- 
stitute 



-i)V 



in (41) to obtain 

CXD 



s=0 



1/2 
s 



r=0 



where 



22"(4n + l)(2r)!(r-f n)! 
{2r + 2n+ l)!(r-n)! 



(42) 



{~iya,r (43) 



P2n{q) (44) 



and we made the substitution 



^2r 



E 

n=0 



22"(4n-f l)(2r)!(r + n)! 
{2r + 2n + iy.{r - ny. 



P2n{q). 



By exchanging the summations over r and n and then n 
and s it is possible to rewrite Eq.(^3|) as 



r{e) - E C2nP2n{0), 



(45) 



C.„ = 22"(4n+1)E( y^)(-l) 
s\ (-l)'^(2r)!(r + n)! 



E 



r J (2r-f 2n+l)!(r-n)!' 



(46) 



FIG. (18) gives the coefficients C2n for n = 1,...,10, 
and for several values of v. While FIG.(|l8|) shows the ex- 
ponential convergence of algorithm with n, it also shows 
that the coefficient of the convergence is small for v ^ 0.9. 
It can be seen that the number of required terms ap- 
proaches 20 for V = 0.9 if the error is required to be less 
than 10~^. (The general sum would have polynomials 
of odd as well as even order, and for each 1, a set of az- 
imuthal quantum numbers spanning -2n to 2n). Hence 
in general, to compute the distorted apparent horizon 
would take a search over 20^ parameters in a minimiza- 
tion routine. This is equivalent to inverting a full ma- 
trix of this size, and would be expected to be slow. The 
boundary value problem is expected to be much faster. 
It is a fact that the boundary value problem as now im- 
plemented does not handle multiple black- hole spaces, so 
its speed is counteracted by the impossibility of using it 
in 2-hole cases. However, it may be possible to use a flow 
method which does recognize the existence of seperate 
black holes, run down to find the two holes each with 
some accuracy, and then to use this boundary soluion 
code to quickly get a highly accurate result. We are con- 
fident such a combined tool would be of great utility. 
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VII. APPENDIX A: EVALUATION OF THE 
RESIDUAL 

The evaluation of Cartesian derivatives on S is car- 
ried out by constructing 3D finite difference stencils at 
each mesh point on S. The finite difference stencil, de- 
noted by M, consists of 26 additional points around each 
mesh point. These 26 points are ±Sx, ±Sy and ±Sz away 
from the central mesh point as shown in FIG.(p^. These 
points, as shown, are organized into three planes of con- 
stant z: z = zq — Sz, zq, zq + Sz. Each plane contains the 
nine nearest neighbors to the center point, including the 
center point itself in the case oi z = zq. We use a single 
discretization scale h {dx = 5y = 5z ~ h) which is always 
proportional to the mesh spacing 56 = tt/{Ns — 1). 

To define tp{x,y,z) at each stencil point x e A/" we 
use its split into radial and angular parts, (p{x,y,z) — 
r — p{9, (j)). For each stencil point x we compute the cor- 
responding spherical coordinates (r,j;,6x,4>x)- This point 
can be thought of as a ray emanating from the origin of 
our spherical coordinate system (which coincides with the 
origin of our Cartesian coordinate system) along [Ox, (px) 
of length r^. FIG.(^o|) labels the point x as P. The dashed 
line from P to the origin is the ray from the origin. Its in- 
tersection with S is denoted by a filled square. The value 
of at X can be obtained by computing (px) via bi- 
quartic interpolation where the truncation error has a 
leading order term which is fourth order in the grid spac- 
ing h. The interpolation is carried out with values of 
p defined on mesh points of S using a 16 point stencil. 
FIG.(pTl) shows the choice of these stencil points in the 
interior of the mesh. At the poles a special choice is made 
of stencil points which takes into account the indentifi- 
cations made at the poles. FIG.(p2|) shows a choice of 
stencil points for an interpolation point near the pole. 
This approach leads to a fourth order truncation error in 
p{Ox, 4>x) at all points on S. Then ip can be constructed 
for every x£Afasip = rx — p{Ox, 4>x)- Using this ap- 
proach (p is defined at any finite difference stencil point 
for every mesh point on S. The finite difference expres- 
sions for A^iy9, A'lA'jip (corresponding to first and second 
derivatives) are computed at each of the mesh points. 
The residual is then evaluated on S using these finite 
difference approximations for the derivatives to 0{h^), 
and metric data {"^ij, dkjij, Kij) which are specified ei- 
ther analytically or interpolated from an enveloping 3D 
Cartesian grid. 

Because we use 0{h?) finite difference approximations 
A^(p, A'lAjip to the derivatives, this approach leads to 

an 0{h?) truncation error in evaluating F[p\ . Because 
of our special attention to points near the pole, F is eval- 
uated smoothly everywhere on S. 

With a means for evaluating F at any point in the 
domain of S it is straightforward to generate J^p numer- 
ically using Eq. (^ . The algorithm for this is summarized 
as follows: 

Specify metric data everywhere on S 
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Evaluate F[p] everywhere on S 
For each point(labelled by D) in S 

Perturb pi, = pu + e 

Specify metric data on perturbed point 

Evaluate + e] at the /2-th point. 

Compute the pv component of the Jacobian matrix 

using (|) 
End loop over points on S. 

This gives the Jacobian matrix, J^p, for F evaluated. 
J^p is a (A^2 _ + 2) x (iV^ _ + 2) matrix which 
is used in Newton's method as foUows: 

Start with an initial guess surface p = po 
while ||F|| > stopping criterion 

Compute the Jacobian J^p for the current p 

Evaluate F[p\ 

Solve J ■ Sp = —F[p] for Sp 
Update the surface p = p + Sp 
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FIG. 2. This figure shows the sparsity structure of the Jacobian matrix for a 33 x 33 run. 
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FIG. 3. This graph shows the log of the error norm of the numerical solution, ||p — p\\2 versus log{Ns). The slope of the 
graph is —2.11 and hence the solution is 0{h?). 
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1 1.5 2 2.5 

Mask Radius 

FIG. 5. Graph of percentage relative error in p versus mask radius. Note that past Rm, = 2.5 the solver did not meet the 
stopping criterion of 10~*. 
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FIG. 6. Percentage relative error as a function of the offset 5. 



18 




FIG. 7. The following figure shows the cross-section of apparent horizons located for various values of a. The solid line shows 
the = TT slice of the apparent horizon. The dashed line shows the 9 = it/2 slice of the horizon. As expected the </> = tt slices 
show increasing deformation for increasing a. 
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FIG. 8. ||f - r+||2 versus A^^ for a = 0.9. 
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Angular Momentum parameter a (J/M) 



FIG. 9. The percentage error in the area of the apparent horizon for a Kerr hole versus a is shown for Ns = 17, 25, 33, 41, 49, 65. 
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FIG. 10. This figure shows the y-z cross-section of Schwarzschild black hole apparent horizons located for boost velocities 
of u = 0.0, 0.1, 0.9 in the y-z direction. The dashed circle in each of these figures is the apparent horizon for an unboosted 
Schwarzschild black hole. Note that the point of contact is along the y-z direction orthogonal to the direction of the boost. 
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FIG. 11. This figure shows the error in the areas of apparent horizons found for a Schwarzschild black hole boosted in the 
y-z direction. 
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FIG. 12. This figure shows the error in the radial coordinate location r = 2M for a Schwarzschild black hole apparent horizon 
for V = 0.0, .., 0.9. The black hole is boosted in the xyz-diroction. 

FIG. 13. This figure shows plots of the apparent horizon for the a = 0.0 runs for v = 0.0, 0.3, 0.6 and 0.9 for boosts in the 
xyz-direction. The mesh used to find the apparent horizon is shown from a top perspective. As the boost velocity is increased 
we see that the surface is contracted in the xyz-direction. 
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FIG. 14. This figure sliows the y-z cross-section of Kerr (a=0.9) black hole apparent horizons located for boost velocities 
of t; = 0.0, 0.1, 0.9 in the y-z direction. The dashed circle in each of these figures is the apparent horizon for an unboosted 
Schwarzschild black hole. Note that the point of contact is along the y-z direction orthogonal to the direction of the boost. 
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FIG. 15. This figure shows the error in the areas of apparent horizons found for a Kerr black hole boosted in the y-z direction 
with a = 0.9. 
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FIG. 16. This figure shows the error in the spheriodal radial coordinate location r = r+ for a Kerr black hole apparent 
horizon (a = 0.9) for v — 0.0, .., 0.9. The black hole is boosted in the xyz-direction. 

FIG. 17. This figure shows plots of the apparent horizon for the a = 0.9 runs for v = 0.0, 0.3, 0.6 and 0.9 for boosts in the 
xyz-direction. The mesh used to find the apparent horizon is shown from a top perspective. As the boost velocity is increased 
we see that the surface is contracted in the xyz-direction. 
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FIG. 19. This figure shows 3 pianos of the finite difference molecule M ■ Plane 2 is the center plane with z = zq. The center 
point is denoted by a filled square which is the mesh point {xq, yo, zq). Each of the other stencil points in plane 2 are a distance 
±6x and ±Sy away and are denoted by circles. Plane 1 is aX z = zo — 5z and plane 3 is at « = 00 + Sz. 




FIG. 20. P is the point at which we wish to estimate ip. The dashed line is the radial line from the origin of our spherical 
coordinate system to P. The filled square on the surface is the interpolation point where we evaluate p{Ox, ipx)- 



FIG. 21. This figure shows the choice of stencil points that we use for biquartic interpolation for interpolation points that 
are on the interior of the grid. The interpolation point is labelled by a filled circle and the mesh points that are used as an 
interpolation stencil are denoted by filled squares. 
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FIG. 22. This figure shows the choice of stencil points for biquartic interpolation for interpolation points that are near the 
poles. We view the pole points from a SDimensional perspective where the identifications of points in the ^Direction is taken 
into account. This leads to the special choice of interpolation stencil points as shown. This leads to fourth order truncation 
error in the interpolant at the poles. The interpolation point is labelled by a filled circle and the mesh points that axe used as 
an interpolation stencil are denoted by filled squares. 
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This figure "figl3.gif" is available in "gif" format from: 



http://arXiv.org/ps/gr-qc/0002076vl 



This figure "figl7.gif" is available in "gif" format from: 



http://arXiv.org/ps/gr-qc/0002076vl 



